Notes & questions

postgibbs thinning: “If you want to keep all the samples you have drawn, you have to put the same number as you put before. Or, you can input a multiple of the original number. For example, in this case, you can input 10, 20, 30 and so on, because the original interval was 10. If you type an inappropriate number, the program will stop with a suspicious message.”
The independent chain size corresponds to the number of independent samples that can be seen as independent. If the independent chain size is 3, the statistics (i.e. posterior mean and SD) are equivalently calculated using only 3 independent samples, and obviously, this is too small. This indicates whether you need more samples or not.
Effective sample size: “At least > 10 is recommended. > 30 may be better.”
You can evaluate a sufficient interval for samples to be saved using Independent chain size and Autocorrelations. Two adjacent samples are usually highly correlated because the next sample is drawn based on the current one. When the correlation is still high between distant samples, the dependency-level is also high and the absolute values of the samples should be very similar.
When the distribution of the samples is not normal, there is a difference between “Mode” and “Mean”

Geweke (1992) convergence diagnostic “is a time-series approach that compares the mean and variance of segments from the beginning and end of a single chain.”

\[z = \frac{\theta_a-\theta_b}{\sqrt{Var(\theta_a)+Var(\theta_b)}} \]

“where a is the early interval and b the late interval. If the z-scores (theoretically distributed as standard normal variates) of these two segments are similar, it can provide evidence for convergence… If the chain has converged, the majority of points should fall within 2 standard deviations of zero… The occurrence of the scores well within 2 standard deviations of zero does not indicate lack of convergence, while deviations exceeding 2 standard deviations suggest that additional samples are required to achieve convergence”

What portions of chain is postgibbs using? Last sample and what?

postgibbs_samples file: "The meaning of each column is:
- Column 1: the sequential number for saved samples.
- Column 2: the actual round in which each sample was drawn.
- Column 3: the number of parameters.
- Column 4 or later: actual samples ordered by the position index."

find data/derived_data/gibbs_varcomp -type f -name 'last_solutions' | sed -r 's|/[^/]+$||' |sort |uniq

nohup echo -e 'renf90.par \n 818000 0 \n 20' | /usr/local/bin/thrgibbs1f90 > gibbs.iter5.3v3alt.continue.out & exit

Sampling parameters:

1,000,000 total samples
Burn-in of 50,000 samples removed in postgibbs
Kept every 20th sample during sampling, thinned to every 100th sample in postgibbs

Setup

gibbs_samples <-
  purrr::map2_dfr(.x = rep(c(1:10), 
                           times = 7),
                  .y = rep(c("1", "2", "3alt", "5", "7", "8", "9"),
                           times = 10),
                  ~ read_gibbs_samples(iteration = .x,
                                 region = .y)) %>% 
  select(iter, dataset, everything(), -X1, -X3) %>% 
  purrr::set_names("iter", "dataset", "round", c(1:14)) %>% 
  tidyr::pivot_longer(cols = c(`1`:`14`),
                      names_to = "param") %>% 
  mutate(param = case_when(param == "1" ~ "dir1dir1",
                           param == "2" ~ "dir1dir2",
                           param == "3" ~ "dir1mat1",
                           param == "4" ~ "dir1mat2",
                           param == "5" ~ "dir2dir2",
                           param == "6" ~ "dir2mat1",
                           param == "7" ~ "dir2mat2",
                           param == "8" ~ "mat1mat1",
                           param == "9" ~ "mat1mat2",
                           param == "10" ~ "mat2mat2",
                           param == "11" ~ "mpe1mpe1",
                           param == "12" ~ "mpe2mpe2",
                           param == "13" ~ "res1res1",
                           param == "14" ~ "res2res2"),
         iter = as.character(iter))

gibbs_mce <-
  purrr::map2_dfr(.x = rep(c(1:10), 
                           times = 7),
                  .y = rep(c("1", "2", "3alt", "5", "7", "8", "9"),
                           times = 10),
                  ~ read_gibbs_mce(iteration = .x,
                                   region = .y)) %>%
  left_join(gibbs_samples %>% 
              mutate(iter = as.integer(iter)) %>% 
              group_by(iter, dataset, param) %>% 
              tally(name = "n_samples")) %>% 
  select(iter, dataset, param, n_samples, everything())

gibbs_psd <-
  purrr::map2_dfr(.x = rep(c(1:10), 
                           times = 7),
                  .y = rep(c("1", "2", "3alt", "5", "7", "8", "9"),
                           times = 10),
                  ~ read_gibbs_psd(iteration = .x,
                                 region = .y)) %>%
  left_join(gibbs_samples %>% 
              mutate(iter = as.integer(iter)) %>% 
              group_by(iter, dataset, param) %>% 
              tally(name = "n_samples")) %>% 
  select(iter, dataset, param, n_samples, everything())

Sanity check

Make sure every parameter has the correct number of iterations

gibbs_samples %>% 
  group_by(iter, dataset, param) %>% 
  tally(sort = TRUE) %>% 
  ungroup() %>% 
  distinct(iter, dataset, n)

gibbs_samples %>% 
  distinct(iter, dataset) %>% 
  arrange(dataset, iter)

Monte Carlo error by time series

Parameters with effective sample size lower than 10

gibbs_mce %>% 
  filter(10 > effective_sample_size) %>% 
  DT::datatable(rownames = FALSE)

Posterior standard deviation

Parameters with fewer than 10 independent batches

gibbs_psd %>% 
  filter(10 > independent_batches) %>% 
  DT::datatable(rownames = FALSE)

Plot post-burnin and post-thinning samples

plot_gibbs_iter <-
  function(df) {
    df %>% 
      filter(!stringr::str_detect(param, "res|mpe")) %>% 
      ggplot(aes(x = round,
                 y = value,
                 color = param)) +
      geom_line() +
      ggsci::scale_color_ucscgb() +
      theme_classic() +
      facet_wrap(~ dataset, nrow = 3)
  }

plot_gibbs_region <- 
  function(df) {
    df %>% 
      filter(!stringr::str_detect(param, "res|mpe")) %>% 
      mutate(iter = glue("Iteration {iter}")) %>% 
      ggplot(aes(x = round, 
                 y = value,
                 color = param)) +
      geom_line() +
      theme_classic() +
      ggsci::scale_color_ucscgb() +
      labs(x = "Round", 
           y = "Value",
           color = "Parameter") +
      facet_wrap(~ iter, ncol = 2)
  }

Desert

gibbs_samples %>% 
  filter(dataset == "3v1") %>% 
  plot_gibbs_region()

Southeast

gibbs_samples %>% 
  filter(dataset == "3v2") %>% 
  plot_gibbs_region()

High Plains, control

gibbs_samples %>% 
  filter(dataset == "3v3alt") %>% 
  plot_gibbs_region()

Arid Prairie

gibbs_samples %>% 
  filter(dataset == "3v5") %>% 
  plot_gibbs_region()

Forested Mountains

gibbs_samples %>% 
  filter(dataset == "3v7") %>% 
  plot_gibbs_region()

Fescue Belt

gibbs_samples %>% 
  filter(dataset == "3v8") %>% 
  plot_gibbs_region()

Upper Midwest & Northeast

gibbs_samples %>% 
  filter(dataset == "3v9") %>% 
  plot_gibbs_region()

Post-Gibbs sampling summary

Harly Durbin